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ABSTRACT 

We present a new symplectic integrator designed for collisional gravitational A-body 
problems which makes use of Kepler solvers. The integrator is also reversible and con¬ 
serves 9 integrals of motion of the A-body problem to machine precision. The integra¬ 
tor is second order, but the order can easily be increased by the method of lYoshidal 
We use fixed time step in all tests studied in this paper to ensure preservation of 
symplecticity. We study small A collisional problems and perform comparisons with 
typically used integrators. In particular, we find comparable or better performance 
when compared to the 4th order Hermite method and much better performance than 
adaptive time step symplectic integrators introduced previously. We find better perfor¬ 
mance compared to SAKURA, a non-symplectic, non-time-reversiblc integrator based 
on a different two-body decomposition of the A-body problem. The integrator is a 
promising tool in collisional gravitational dynamics. 

Key words: celestial mechanics- gravitation- methods: numerical- methods: 

analytical- globular clusters: general. 


1 INTRODUCTION 

Astrophysical Abody problems can be classified by their 
collisional nature. We define collisionle ss problems as those 
where the two body relaxation time fsee lBinnev fc Tremaine! 
(hoosB is larger than the timescale of interest. A typical 
globular cluster has t re iax ~ 400 Myrs, whereas a dark mat¬ 
ter (DM) halo’s relaxation time is many orders of magnitude 
larger than a Hubble time. If we are concerned with evolu¬ 
tion over the time scale of the Universe, the former is colli¬ 
sional while the later is collisionless. In the collisionless case, 
a test particle’s motion is approximately due to a smooth 
potential of the mass distribution. In this case the dynamics 
is given by the collisionless Boltzmann equation, a partial 
differential equation in 6 variables. Collisional problems are 
much more complicated, requiring a solution of the 6A di¬ 
mensional Liouville equation, with A the particle number. 

DM is typically believed to be collisionless. Collisionless 
problems, like those in DM simulations, are typically solved 
by a Monte Carlo sampling of A particles of the phase fluid. 
To avoid non-physical two-body encounters and two-body 
relaxation, researchers usually include a softening length 
such that the force goes to a finite con stant as i nter par- 
ticle distances go to 0 (see for example IVogelsberger et all 
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d2008lU . The dynamics then remain governed approximately 
by a collisionless Boltzmann equation. H owever, to o high a 
softening length leads to other errors dBertschingerj[l998l l. A 
study of errors in dynamics int roduced by softening lengths 
is carried out in lDehnenl d200ll l. The impact of the softening 
length cannot be overstated: without this feature, modern 
codes would not be able to carry out the large DM A-body 
integration with A > 10 9 . 

We cannot include a softening length in collisional prob¬ 
lems, such as in globular clusters, where close particle en¬ 
counters drive the dynamics. Tight binaries, hierarchical 
triples, and higher order systems will form. In order to model 
these systems accurately, standard integrators will require 
computationally prohi bitive resources in the form of small 
time steps dKvaerno fe Leimkuhled l200dl . These new sys¬ 
tems also introduce new dynamic timescales, and standard 
integrators will not be suited for the resulting large dynamic 
range of timescales of the problem. 

The most popular work aroun d t hese p roblems is regu¬ 
larization, described in lAarseth et all (|2008h . The basic idea 
is to perform a coordinate and time transformation such that 
in the new coordinates, the equations of motion remain non¬ 
singular as the inter-particle separation goes to 0. There are 
ways to regularize three (and more) body encounters. 

With or without regularization, most integrators up un¬ 
til now show inaccuracies when applied for long times to 
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general collisional systems. One way of mitigating_this prob¬ 
lem i s to incorporate a geometric al integrator (lHairer et al.l 
120061 : ILeimkuhler fe Reicbl 120041 ') , an integrator that con¬ 
serves qualitative features of the original problem. Such fea¬ 
tures can be due, for example, to continuous symmetries 
in the Hamiltonian, leading to integrals of motion, or dis¬ 
crete symmetries such as time reversibility. Symplectic in¬ 
tegrators are one such geometric integrator which preserve 
phase space volume. The advantages of such an integrator 
can be powerful for long term evolution of problems due 
to a lack of secular drift in energy errors and preservation 
of global structure in the governing differential equations. 
For this reason t hey are the pre ferred method in cosmolog- 
ical simulations ISprjnge] 120051 ) and planetary simulations 
dWisdom fe Holmanl 1991 ), but have not been successfully 
employed yet in globular clusters or other collisional sys¬ 
tems. 

Symplectic integrators have not found success in their 
application to collisional systems in the past because when 
time steps are adapted as a function of phase space or if 
the integrator changes form as a function of pha s e spa ce 
(described as a “switch” in iLeimkuhler fe Reichl (l2004r) f. 
the integrators no longer approximately obey a Hamilto¬ 
nian, a major feature which hel ps make the m success ful in 
the first place. Non e theles s, iMikkola fe T anikawa (199^) and 
IPreto fe Tremainel (1999) independently discovered a sym¬ 
plectic integrator for small N collisional problems which is 
able to adapt its time step as a function of the total potential 
energy of the system. 

Time reversible integrators share similarities to some 
symp lectic integrators for regular or near regular mo¬ 
tion dHairer et all 120061 ) and indeed time reversible meth¬ 
ods for collisional systems have been developed. Espe- 
cially notewor t hy is a n adaptive symmetric Hermite method 
iKokubo et al.| (1 19981) impl emented in the popular Aarseth 
codes ( Aarseth et all 120081 1. but the computational cost of 
the integrator is high and its use is recommended only for 
planetary systems. The leapfrog method has also been made 
to be adaptive and t ime reversible: this integrator is studied 
in iGoncalves Ferrari et alj d2014T) and other works, but its 
performance is not found to be as favorable as other existing 
methods. 

In this paper, we present a new symplectic integrator 
for collisional gravitational 7V-body dynamics. The integra¬ 
tor is inspired by the non-symplecti c and non-reversible in¬ 
tegrator in iGoncalves Ferrari et ahl d2014l b SAKURA, and 
makes use of Kepler solvers. Like SAKURA we decompose 
the TV-body problem into two-body problems. In contrast 
to SAKURA, our two-body problems are not independent. 
The integrator is reversible and symplectic and conserves 
9 integrals of motion of the IV-body problem to machine 
precision . The integrator is second order, in the sense of 
lHairer e t al.| (12006 ) . but the order can be increased by the 
method of lYoshidal (119901) . We use fixed time step in all tests 
studied in this paper in order to remain exactly symplectic. 
We study small N collisional problems and perform compar¬ 
isons with typically used integrators. In particular, we find 
comparable or better performance when compared to the 
4th order Hermite method and much better performance 
than the second order sym plecti c integrator introduced b y 
IMikkola fe Tanikawal (1 19991 ') and IPreto fe Tremainel lll999lb 
We also find improved performance over SAKURA. Thus 


our integrator is promising as a tool in collisional gravita¬ 
tional dynamics. We plan larger N tests of the method in 
future work. 

The organization of the paper is as follows. In Section 
O we outline the method along with its general properties 
and give first illustrations of its power in handling collisional 
problems. Section 12.21 provides a mathematical derivation of 
the method along with a proof of its symplecticity and other 
properties. Section [2.2.11 shows numerical evidence for sym¬ 
plecticity and the existence of a Hamiltonian closely followed 
by the method. 

Section [3] is dedicated to numerical tests. Section 0 
studies tests of regular 3 body systems in order to clearly 
show the advantages of the symplecticity of the method over 
other commonly used integrators. In Section 13.21 we study 
the more common chaotic IV-body problem and show the 
good performance of the method persists over other integra¬ 
tors. We also argue why symplectic methods are still highly 
desirable for chaotic problems. Finally, Section []] shows how 
to combine the method with a traditional leapfrog method 
while maintaining symplecticity and reversibility. 


2 THE METHOD 

2.1 Basic use and overview 

First we outline the steps for the basic integrator in pseu¬ 
docode, a map of phase space coordinates <f >It is a sym¬ 
plectic, symmetric ( or reversible) second order integrator, 
although i t is easy to ge neralize to other orders via the 
method of lYoshidal d 19901) . First we begin with a first or¬ 
der mapping, <f>h■ A drift is defined for a coordinate system 
in which the A-body problem is separable: the kinetic en¬ 
ergy depends only on momentum and the potential energy 
only on coordinates. For example, a drift in cartesian coor¬ 
dinates is x' = x + hv, where x' is a new position, and x 
and v are the old position and velocity respectively, and h 
is a time step. 

Drift all particles for time h ; 
for all pairs of particles (i, j) do 
Drift particles i and j for time —h ; 

Apply a Kepler solver to advance the relative 
coordinates of i and j by h ; 

Advance center of mass coordinates of i and j by 
h ■ 
end 

Algorithm 1: A Kepler integrator mapping f>h- 


Now we define an integrator that has the same steps 
but exactly in reverse as (jr h . Our second order integrator is 
then their composition: 

<Ph = <t’h/2<f>h/2 (!) 

There is no preferred order of particles in the loop of Algo¬ 
rithm [T] in general. The order probably affects the specific 
form of higher order terms in the Hamiltonian obeyed by 
the method (described in Section 12.21) we did not calculate, 
but we did not find significant differences in accuracy when 
changing the order in which particles are updated for dif¬ 
ferent problems. There could be special problems for which 
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a specific ordering gives better results. The composition of 
the two methods makes the map time reversible, as will 
be discussed below. With the two compositions there are 
A(A — 1) total Kepler solvers, where A is the number of 
particles. These Kepler solvers account for the only signif¬ 
icant computation time. It is possible to create a second 
order integrator with half as many Kepler solvers (in fact, 
Algorithm |T] is such an example), but this is the simplest 
we could find that is symmetric, and thus time reversible. 
It can be shown theoretically that a reversible integrator 
alone already shares many properties of a symplectic inte¬ 
grator when applied _to integra ble and near-integrable sys¬ 
tems (Chapter Xl lHairer et al.l J200(j )'). While the A-body 
problem is not integrable or near-integrable, we can expect 
advantages in structure preservation by using a reversible 
method. 

We can sho w sim ply <jf i s rev ersible by checking def¬ 
inition V.1.4 in lHairer et alj d2006l) . For a reversible map 

<t>h, 


(fthfi—h — F, 


( 2 ) 


where I is the identity matrix. A numerical method should 
obey eq. [2] independent of the step size h, to round off er- 
ror(j We find has such a property. We note here that, 
contrary to their claim, we do not find the method of 
iGoncalves Ferrari et al.l ( 20141 '). which also provides an inte¬ 
grator with Kepler solvers, to be time reversible. Tests show 
that the error in eq. [2] applied to their method grows with 
h and is not at the level of roundoff. We prove SAKURA’s 
irreversibility and explain why SAKURA is not reversible in 
Appendix m 

When compared to a popular Drift -Kic k- Drift (DKD) 
l eapfr og method (see the introduction of lPreto fc Tremaine! 
(4999 ) for a description of DKD and its counterpart KDK), 
as far as the main computational work is concerned, we have 
substituted force calculations between the A(A — l)/2 pairs 
for solving A(A — l)/2 Kepler problems. Generally, the Ke¬ 
pler solvers will take slightly longer than the force calcula¬ 
tions. We can compare the computational cost of solving 
a two body problem (2BP) between 4> h and DKD. For the 
2BP, <ff reduces to a two body Kepler solver, assuming the 
center of mass is stationary. A two body Kepler solver ad¬ 
vances 6-dimensional relative phase space coordinates. 

To test this, we choose reduced mass fi = 1, total mass 
M — 1, eccentricity e = 0.9, semi-major axis a = 1 (so that 
the period is 27 t), and time step h = 0.01 in units where 
Gil. The initial conditions are at apoapsis: x = 1 + e, 
y = 0, v x = 0, v v = y/(l — e)/(l + e). Running the 2BP 
for 25 orbits on a Macintosh laptop computer, the compu¬ 
tation time for this problem is about 38% longer for the Ke¬ 
pler solver (for equal h) compared with DKD. Our Kepler 
solver uses universal variables for advancing the gene ral Ke¬ 
pler p roblem wi th bound or unbound orbit s following iDanbvl 
(Il988h : see also lMikkola fc Innanenl (Il999fl . 

Map (j>\ is symplectic as we will prove analytically and 
show numerically. We will also show numerically that of the 
10 integrals of motion known for the general A > 3 A-body 


1 Actually, for several reversible integrators we tested, the error 
grows for h > 1. This is due to errors in the finite precision 
arithmetic and a loss of decimals stored. 


problem, conserves 9 to machine precision. As for the last 
integral of motion, the energy, a surrogate Hamiltonian H 
is conserved that differs from the A-body Hamiltonian by 
the order h 2 , because it is a second order integrator. It will 
remain close to the energy as long as h is small enough, es¬ 
pecially in a ne ar regular sys t em. F or rig orous results on thi s 
statement, see lHairer et al.l (l2006l ~) and lEngle et al.l ((2005) . 
The DKD leapfrog integrator, in addition to being symplec¬ 
tic and reversible, also has the same conservation properties. 
However, it is not as well suited for collisional studies as 
the explanation for this lies in the behavior of the infinite 
series in h found in its H\ specificall y, theorems for long 
term energy conservation presented in lHairer et al.1 (|2006n 
no longer apply. 

We illustrate this point in Table |T| Here we test the a 
three body problem with the two integrators, DKD leapfrog 
and (pfi- We choose the three body Pythagorean problem 
(3BPP), described in detail in ISzebehelv fc Peters! (|l967l l. 
with their same initial conditions. We use time step h = 
0.01. For reference, by time t = 10, there have been 6 close 
encounters. The first close encounter happens shortly after 
t = 1.5. 

We run the integrators for t = 2.0, a time shortly af¬ 
ter the first close encounter. Even when (j>\ is given a time 
step 15 times larger, and its computation time is more than 
5 times smaller, it performs similarly to a DKD integra¬ 
tor. This test shows the unsuitability of DKD for collisional 
studies. We see conservation of all integrals of motion of this 
problem: the energy, angular momentum, and four center of 
mass coordinates. Section 12. 2 1 explains the details of why (jf 
has these conservation properties. We will revisit in Section 
I3.2l the unsuitability of leapfrog for use in chaotic problems. 


2.2 Derivation of method and basic properties 

This Section and Section 12.2.11 derive and prove the conser¬ 
vation properties of Readers not interested in the tech¬ 
nical details may wish to skip to the performance of the 
integrator when compared to other integrators in Section [3] 
To derive (ff , first proceed as in Section 12.11 We first 
build a second order integrator cj>h from Algorithm |T] We 
will follow notation that is conventional in lYoshidal ll 199(f) . 
lHairer et al.l (120061) . and others. We motivate the integrator 
by writing first the A-body Hamiltonian, 

H = T + V 

T + (3) 

= r +EE(^- T ' 3 )- 

i j>i 

Here T is the kinetic energy, V is the potential energy from 
Newton’s gravitational law, Vij is the pairwise potential of 
particles, Tij is the kinetic energy of particle i plus particle 
j, and Kij is the two body Hamiltonian for particles i and 
j (K stands for Kepler) including both relative and center 

of mass terms. _ 

Using an operator splitting method (see lYoshidal 
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Table 1. Comparison of performance of two second order methods: DKD leapfrog integrator vs for a collisional study of the 3BPP 
described in the text, tmax is 2.0, shortly after the first close encounter. Even when is given a time step nearly 10 times larger, and 
its computation time is more than 5 times smaller, it performs at least as well as a DKD integrator. 



\AE/E 0 \ 

dL 

dx cm, 1 

^cm,2 

dp cm, 1 

dp cm,2 

^cpu 

h 

DKD leapfrog 

8.2 x 10 -6 

1.2 x lO” 13 

8.4 x 10“ 15 

1.3 x 10- 13 

2.0 x 10- 13 

2.1 x 10“ 14 

7.5 

0.0001 

<t>l 

3.7 x lO" 6 

1.1 x 10 -13 

4.7 x 10 -14 

1.4 x 10 -14 

2.6 x 10 -14 

8.0 x 10“ 15 

1.4 

0.0015 


GUI)) , we construct 


<ph=\ fl exp (hD Kij ) exp (-hD Tij ) I exp (hD r ) . 

\(i,j) pairs 1 

n 

( i,j ) pairs J 

(4) 

where each ip mapping is defined as one of the exponentials. 
Here Da is defined by Da} = {/, A}, where {} are Pois¬ 
son brackets. Exponentials are calculated by writing their 
Taylor series. Thus we can check that mapping iph corre¬ 
sponds to a drift of time h, defined in Section mu Equation 
[4] corresponds exactly to Algorithm [l] As mentioned before, 
we found no generally preferred ordering in the product of 
eq. ©• Note that the ip in parenthesis are invariant under 
exchange of i and j. 

First we prove the symplecticity of map ©■ Because a 
composition of symplectic mappings is symplectic, we will 
show each exponential mapping in eq. © is symplectic. In 
the canonical coordinate basis discussed below, the symplec¬ 
tic matrix is defined by 



ft = 


0 

-I 


I 

0 


I and 0 are the 3N x 3 N identity and zero maps respectively. 
For a map with associated Jacobian J and Jacob i an trans¬ 
pose J T , the symplectic condition says (ISussman fe Wisdoml 
1200111 


ft = J T fiJ. 


(5) 


We check this condition for eq. ©. starting with ip We 
can calculate the Jacobian resulting from this map acting 
on a phase space y, where y is a 6 N vector composed of the 
coordinates and momenta, in that order. We use cartesian 
coordinates and their associated linear momenta from the 
Wbody Lagrangian. J is given by 


I M 
0 I 


( 6 ) 


Here, M has the property 

M T = M. (7) 


We can now verify the symplectic condition ©: 


J T SiJ 



I 

M T -M 


This proves the familiar result that a drift of all particles is 
a symplectic map. 


Next we focus on ip ^ 3 , a drift of particles i and j. The 
associated Jacobian again has form eq. © with property © 
so map 1 pf^ 3 , a drift of particles i and j is symplectic. 

Finally, we have the map ipfp 3 . By definition, this map¬ 
ping is a Hamiltonian flow of particles i and j and thus sym¬ 
plectic with respect to their 12 dimensional phase space. It 
is easy to show then using condition © that the Jacobian 
of the entire phase space vector y is also symplectic. Since 
we have shown ipip^ 3 , and ipf^ 3 are symplectic maps, we 
have proven (ph is a symplectic map. 

We discuss the integrals of motion associated with the 
continuous symmetries of the Hamiltonian obeyed by c ph- 
We are first able to show that the total angular momentum 
is conserved exactly to machine precision. To do this, we 
can check that ip\, ip ^ 3 , and 1 p^ 3 conserve the total angular 
momentum. The 6 center of mass integrals are also conserved 
exactly by these maps. It is easy to check the conservation 
properties analytically, and numerical experiments support 
the finding. 

DKD leapfrog shares the property with (ph that 9 of the 
10 integrals of motion are conserved exactly, and it is sym¬ 
plectic. DKD is a well known integrator, but let us discuss 
its properties carefully as they will prove useful for Section 
2]. We have already studied the drift maps D and all that 
remains for analysis is the kick map K. 

We can represent K as exp(h.Dy). Dy is a sum of pair¬ 
wise potential energy operators that commute because the 
potential energies depend only on position. Thus, exp (hDy) 
can be written as a product of maps of form exp {hDy ij ). Let 
us study the map exp (hDyD. This map shifts the momen¬ 
tum of particles i and j by a term proportional to their 
mutual force. The Jacobian of K has a form similar to those 
we have seen: 


J = 


I 

M 


0 

J ’ 


(8) 


with M a symmetric matrix. Then, by similar arguments 
to those following equation © , K is symplectic. The con¬ 
servation of total momentum follows from Newton’s third 
law. Because position coordinates remain unaffected by K, 
conservation of the remaining center of mass constants fol¬ 
lows. The total angular momentum is conserved because the 
change in momenta is a function only of position. 

We see that a discrete symmetry of the problem, time 
reversibility, is not respected by (ph . It is obvious from in¬ 
spection of (ph that eq. © does not hold and one can indeed 
verify this numerically by running the integrator forwards 
and backwards for a time step. We should try to correct 
this feature since, as noted previously, time reversibility un¬ 
der some circumstances provides benefits similar to those of 
symplecticity. 
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When one tests the order of tfih, one finds it to be a sec¬ 
ond order integrator. Where does this come from? It comes 
from a simplification in the final integral of motion of 4>h, 
a surrogate Hamiltonian H. We will investigate where H 
comes from now. To begin, we state that the exact Hamil¬ 
tonian canonical coordinate transformation after time t is 
given by 


y(t) = exp ( tD H ) y( 0). (9) 

Eq. ([9]l follows from the Hamiltonian equations of mo¬ 
tion for a time-independent Hamiltonian. Now, we write 
cj) h = exp {hDfj) (note we are using a H instead of H 
operator). Then we use eq. (|4 )l wit h the Baker-Campbell- 
Hausdorff formula (see lHairer et all 1 2006 11 to obtain a re¬ 
lationship between D^, DK i:j , D Fij , and Dt ■ Next, we use 
the following relation, which can be verified directly: 

[D f , D g ] = T1{g,f}, 

with [] referring to commutation. With this relationship, we 
can solve for H\ 

H = H + h/2 (- {V,T} + +<3( ft2 ) 

\ i j>i / 

= H + 0(h 2 ), 

( 10 ) 

i.e. the order h term vanishes exactly. This unexpected re¬ 
sult explains why the order of <f>h is actually 2 numerically. 
Illumina ting theoreti c al an d nu merical work on H can be 
found in lEngle et al.l (l2005l l and lHairer et al.l d2006l l . Even 
though the infinite series in H generally diverges, they show 
the error in conservation of the first N terms in H, with N 
determined by h, is exponentially suppressed over exponen¬ 
tially long times. Because the 0th order of H is the energy, 
the long term conservation of energy follows. 

4>h is second order, and conserves 9 out of 10 integrals 
of motion exactly. But we can do better and require time 
reversibility, albeit at the cost of twice the number of Kepler 
solvers (two Kepler solvers per pair of particles ). To do this 
we first introduce the adjoint map as defined in lHairer et al . 

<2006h : 

<1 = Cl ( 11 ) 

lHairer et al.l (j2006l ^ shows that adjoint methods can be used 
to increase the order of an integrator, but this will not work 
for rf>h ■ Adjoint methods have other uses: for a reversible 
method we have 


Ch — d 

an equation we can see, by inspection, does not hold for ij>h- 
But the map (jf = <t>' h / 2 < t > h /2 will be symmetric (while still 
incorporating the other conservation properties discussed 
above). (j>\ has another property not shared by <f>h in the 
form of H. Because 4>\ is sy mmetric, i ts H will only have 
even terms in h (for a proof, see lYoshidal l|l99(t )l. specifically: 

H = H + f(y)h 2 +g(y)h 4 + ..., (12) 

where / and g are known functions of phase space. This is 
why symplectic and reversible integrators can only be even 
order. 
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2.2.1 Symplecticity and the perturbed Hamiltonian H 

We first wish to study the properties of H. It is hard to an¬ 
alytically calculate / and g from above but an easier alter¬ 
native exists. Consider an alternate symplectic, symmetric, 
second order map, slightly more complicated than 

Ch = Ct/ 2 Ch/2 

with 

Ch= ( n (i3) 

\ (i,j) pairs J 

The difference with (j)h is that there are twice as many drift 
steps. We write the pseudocode for integrator eq.Q?|]for clar¬ 
ity in Algorithmic] Using a Kepler map sandwiched between 


Drift all particles for time h ; 
for all pairs of particles (i, j) do 

Drift particles i and j for time —h /2 ; 

Apply a Kepler solver to advance the relative 
coordinates of i and j by h ; 

Advance center of mass coordinates of i and j by 
h ; 

Drift particles i and j for time —hj 2 ; 

end 

Algorithm 2: Pseudocode for (jk defined by eq. 1131 


two drift maps in eq uation m reminds us in form to an in- 
tegrator presented in Rein fe Tremain5 (l201lh . equation 17. 
I Rein fe TremairO (1201 if) studied the three-bodv problem in 
Hill’s approximation and solved for the relative coordinates 
of the small masses in a rotating frame. 

The symmetric nature of the terms within parenthesis 
makes computing H easier for Cf than for . The result is 
as follows: 


H = H 2 + 0(h ) 


(14) 


with 


H 2 = H + h 2 ( d {{ T, V} ,V} - ±{{V,T} ,T} 


i j>i x 


} i Vn} {{^»j i Tij} , Tij } 


Notice how H has no dependence on the order of pairs in 
the integrator to second order. First, we can verify our result 
that H should only contain terms even in h. By measuring 
the error in H 2 , we can verify it scales as h 4 as we would 
predict. Next, we can verify that H 2 should be better con¬ 
served than H f or sm all enough h; similar experiments were 
carried out in lEngle et al.l (12005l l. To do this, we use the 
3BPP again for time t = 1.5, before the first collision and 
choose h = 0.01. We indeed find H 2 to be conserved better 
than H, as we show in Figure |TJ 

We note the power law increase in global energy error 
in Cf. The initial slope of both curves in Figure[l]is « 2. The 
higher order RK4 shows a slope of « 2.3 for the same time 
step. The slope in Figure [T] is sensitive to initial conditions. 
In particular, we note that since the 3BPP starts from rest, 
the Kepler solvers in the first step will all solve orbits with 
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Figure 1. Plot showing the improved conservation of H 2 as com¬ 
pared to conservation in energy E, as expected for small enough 
h. We use the map and its associated II 2 , applied to the 3BPP 
run for time t = 1.5 with h = .01, before any collisions occur. 


approximately parabolic motion. If we perturb the veloci¬ 
ties in the initial conditions such the orbits are elliptical or 
hyperbolic, we find the slope of fa to be ss 1 in all cases, i.e. 
the error is proportional to the number of time steps. The 
same holds for the RK4 integrator: the slope is sa 1 for all 
sets of initial velocities we tried. 

Linear growth in global energy error is typical for non- 
symplectic integrators as we will also see in Figure [3] be¬ 
low. The fact that we see similar error growth in fa is 
due to the chaotic natur e of the problem and is expected 
dChannell fc ScovelHl990l '). This is not necessarily a cause for 
concern, and (fa still yields very positive results for chaotic 
problems. We will explain why symplectic integrators can 
still be well-suited for chaotic problems and study chaotic 
problems in Section T3. 2 1 

Now, we return to the simpler q')f l . For all tests we 
tried, fa was less efficient than (fa, so we do not consider 
it further. We can numerically demonstrate the symplectic- 
ity of (fa. To do this, we apply (fa for one time step to the 
3BPP and calculate a J acobian using Richardson extrapo¬ 
lation (|Press et alJl2002h . We then calculate the right side 
of eq. © and subtract it from the left side. This quantity 
should be 0 to roundoff error for a symplectic method. We 
take the absolute value of this difference matrix and sum all 
the array elements, and call this quantity dl. We compare 
the result in dl for (fa, KDI< leapfrog (we find similar re¬ 
sults for DKD leapfrog), a Runge-Kutta 2nd order method 
(RK2), and SAKURA. D stands for drift, defined previ¬ 
ously, and K stands for kick. A kick is defined, as in the case 
of a drift, for a separable Hamiltonian: p' = p + hf, where 
p’ is a new momentum, and p and f are the momentum 
and force respectively. We then calculate dl' = dl/m, with 
m the number of maps per integrator. For example, (fa has 
m = 14. 

The results are shown in Table [2] Even though we cal¬ 
culate dl for only one step, we see a clear difference be¬ 
tween the top three symplectic methods and bottom 2 non- 
symplectic methods. 


Table 2. Symplecticity (measured by dl') comparison of fa with 
other second order methods described in the text after one step in 
the 3BPP. The three known symplectic methods show the small¬ 
est values of dl'. 


dl' 

KDK 

2.0 x 10- 13 

fa 

2.8 x 10~ 13 

4 KDK steps 

7.9 x 10- 13 

SAKURA 

7.4 x 10 -11 

RK2 

1.8 x lO" 9 


3 NUMERICAL TESTS OF METHOD 
3.1 Regular systems 

We now test the performance of (fa against standard in¬ 
tegrators when applied t o A-bo dy pro bl ems. SAKURA is 
presented in Listing 1 of lGoncalves Ferrari et akl (|2014l 'l as 
Python code. To remain faithful to Listing 1, we imple¬ 
ment it and the other codes in Matlab, which has similar 
syntax to Python. We first study regular problems, where 
the properties of the symplectic integrator are most clearly 
seen. The importance of regular problems in testing sym¬ 
plectic integrators is significant. For small h, the KAM tori 
characteristic of Hamiltonian systems translates into close 
new invariant KAM tori u nder symplectic discretization (see 
IChannell fe Scovell rtl99d 'lh The implications of the near in¬ 
variant KAM tori can be seen when we study special peri¬ 
odic solutions of the A-body problem, as excellent long term 
behavior for small enough h. 

In standard collisional codes, even those with regular¬ 
ization, the lack of a symplectic integrator means there will 
be secular error growth in the integrals of motion, linear 
in time for many common methods. The growth is charac¬ 
teristic of the integrator, and distinct from the growth of 
roundoff error. 

For our first problem, we study a hirerchical triple prob¬ 
lem (HTP). In this system there is a tight binary initially 
aligned with the s-axis and at apoapsis. In units where 
G = 1, these particles have mass m 1 = m 2 = 1/2. Their 
semi major axis is ao = 0.01 and eccentricity is eo = 0.9. 
Their orbital period is therefore Pq ~ 0.0063. Their initial 
velocities are oriented in the y direction. 

There is a third particle with m 3 = 1 that forms a bi¬ 
nary pair with the center of mass of the first pair of particles. 
This second pair also begins aligned with the x -axis and at 
apoapsis. The semi-major axis is ai = 1. The orbit is circu¬ 
lar, so the period of this second system is Pi m 0.14, which 
is over 20 times Po- The initial velocity of this pair is also 
in the y direction. The center of mass position and velocity 
of the triple system is 0. The total run time for this system 
will be t = 1.4, which corresponds to about 10 periods of the 
larger binary and 220 periods of the smaller binary. During 
this time, any non-symplectic or non-geometric integrator 
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t/P 

Figure 3. Energy drift for different i ntegrators when applied to 
the f igure-eight three body problem (|Chenciner Sz Montgomery! 
hooch . We choose different time step parameters for the integra¬ 
tors such that the initial energy error is comparable. There is no 
drift for the symplectic methods. 

will show dissipation and secular errors in the integrals of 
motion. 

In Figure [2] we compare the conserved quantities for </>f, 
SAKURA, and a fourth order Hcrmite method. We choose 
time step parameters such that the cpu time of the integra¬ 
tors is approximately the same (t cpu ~ 50). Hermite uses 
n = 1 0 ~ 3 (?7 is a proxy for the time step, see lAarseth et ahl 
1120081) ). <f>\ uses h = 1.5x 10 -5 as does SAKURA. The figure 
shows the result that no secular growth occurs in any inte¬ 
grals of motion for and the error in energy is bounded. 
4>h performs better by all proxies compared to the other in¬ 
tegrators. 

In Figure Owe saw the first signs of the bound on en¬ 
ergy error due to the combined use of reversible and sym¬ 
plectic methods. Another periodic solution is the figure- 
eight three body pro blem , des cribed mathematically in 
IChenciner fc Montgomery! (120001 ). This time our goal is 
to test how important are the properties of symplectic- 
ity and time reversibility in long term behavior. We an- 
alyze this problem with the s ame initial conditions as 
IChenciner fc Montgomery! ll2000l ) and study its evolution 
under different integrators. For this problem, we choose dif¬ 
ferent time steps (and ?? parameter in the Hermite code) for 
the integrators such that the initial energy error is compa¬ 
rable. 

We run the problem for 100 periods. The period of the 
problem is approximately T = 6.32591398, and in this time 
there are 7 particle crossings at the origin. Energy errors are 
averaged over the period to minimize periodic, expected, 
variations in energy conservation when plotting. The re¬ 
sult is shown in Figure [3] The symplectic methods all have 
bounded conservation properties over long periods of time 


as expected, (fa is the non-reversible integrator in Algorithm 
|T]we used to derive cfa. Despite (fa being non-reversible, its 
good behavior remains. SAKURA shows a linear drift in 
energy error, as does the Hermite integrator with a larger 
slope. They, like RK2 or RK4, have an energy error propor¬ 
tional to the number of time steps taken. 

3.2 Chaotic systems 

Next we show the effect of (fa when applied to chaotic sys¬ 
tems, the more common astrophysical scenario. Long term 
conservation of energy is sometimes not apparent in chaotic 
systems when a symplectic integrator is used, as we saw 
in Figure [ 1 ] This is due to the changing nature of the in¬ 
finite series in H. Thus, a naive analysis based on energy 
conservation alone may make one believe that a symplectic 
integrator performs similarly to standard integrators when 
ch aos is involved. Thi s analy sis would be false. As explained 
in lChannell fc Scovell d 19901 ). symplectic integrators in this 
case are still highly desirable and preserve structure of the 
topology of the differential equations by correctly avoiding 
stable invariant objects. We can see this result by applying 
the symplectic Euler map to the simple pendulum problem. 
The map in this case is equivalent to the Standard M ap (see 
a disc ussion and definition of the Standard Map in lYoshidal 
(ll993l) ). While orbits in the Standard Map can be chaotic, 
for small enough h they avoid invariant topological objects 
and the chaos is bounded._ 

Also relevant are results in lMcLachlan fc Atelal (1 1992li . 
where they applied symplectic integrators to chaotic prob¬ 
lems. A particularly interesting result is that they find 
that the long-term statistics of behavior of orbits con¬ 
verges quite rapidly for the perturbed linear oscillator, a 
chaotic problem. This reminds u s of a recent result by 
IPortegies Zwart fc Boekholtl d2014l) where it is found that 
statistics of orbits for three-body problems converge to cor¬ 
rect answers for accurate enough integrators. An interesting 
question is whether the convergence rate improves for a sym¬ 
plectic integrator like (fa , which they did not test. 

We proceed with these results in mind. Binaries in glob¬ 
ular clusters and other collisional astrophysical systems have 
been the subject of much study and pose one of the great ¬ 
est challenges to the se simulations (Hegg ie fe Hutl d2003l) . 
lAarseth et alJ (120081 )). Realistic globular clusters have a 
fraction of stars in primordial binaries and here we study 
this effect again in a toy problem. We choose a problem of 5 
bodies sampled from a Plummer sph ere. The properties o f 
the Plummer sphere are described in [Heggie fc Hull I 2 OO 3 I) . 
Standard units are used in which G = M = 1 and E = 1/4. 
M is the total mass. Each particle has mass m = 1/N. This 
corresponds to a half mass radius Th ~ 4/5. In these units, 
the crossing time is f C rosa ~ 1.1. 

Now, we add an extra body that forms a tight binary 
with one of the original 5 bodies. The binary has semi major 
axis a = 0.01 and eccentricity e = 0.9. The mass of each par¬ 
ticle in the binary is 1/10. In practice this means the binary 
is twice as bound as the Plummer sphere (E = —1/2). We 
run this problem for t/taoss = 4.5 with various integrators 
we have discussed, including t he 4th order ge neralization of 
leapfrog labelled as ‘Yoshida’ dYoshidalll990l ). The result is 
shown in Figure [4] (fa outperforms the other methods in 
conservation properties for a given level of computing effort. 
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Figure 2. The HTP problem, described in the text, run for t = 1.4. The three integrators are normalized to approximately equal cpu 
effort. </>? has the best performance across all integrals of motion and its property that it has a conserved, nearby H, is seen by the 
periodic behavior in energy conservation. 


So fyh. is an integrator which shows good performance when 
primordial binaries are present. For binaries sufficiently tight 
and far away from the other particles, <p\ will solve their mo¬ 
tion more accurately for a given compute time than the other 
integrators. As mentioned in Section 12.11 we can increase 
the order of our integrator from a second order integrator 
by the method of lYoshidal (1 19901 1. For this problem, tests did 
not indicate significant accuracy improvements for a given 
computation time using the fourth order reversible general¬ 
ization of <j>\. For general Plummer sphere problems without 
any binary formation we find comparable performance be- 
twee n (pj and a Fle rmite method. H owever, based on results 
from iGoncalves Ferrari et al.l (l2014l l with SAKURA, where 
they found their method reproduced the cluster radius accu¬ 
rately and more rapidly than other methods for a moderate 
N Plummer model, we expect (p^ to perform well for larger 
Plummer models. This is especially true since (p\ shows im¬ 
proved performance over SAKURA in all our tests. 

We now revisit the 3BPP presented in section [HU The 
purpose here is to compare the perfor mance of <jyf,, with the 
second ord er symplectic integrator o f lMikkola fc Tanikawal 
dl999f l and IPreto fc Tremaina (| 1999h , which we Call (pTiem, 
and leapfrog . The time step funct ion for ^Trem will be that 
presented in lMikkola fc Tanikawal d 1999T) . Section \2 . 2 1 shows 
leapfrog’s conservation of all integrals of motion, except for 
the energy, to machine precision. <j >Trem also conserves these 
integrals well, so we will only focus on energy conservation. 


We consider the conservation in energy for a given computa¬ 
tional effort (the computational effort is changed by varying 
the timestep). We run the 3BPP for t = 4: this corresponds 
to a time shortly after the third encounter. The result is 
shown in Figure [5] 

We see the poorest integrator is DKD leapfrog, which 
is not surprising. Improved performance is shown by cp Trem, 
but the best results are obtained by (p\. The superior per¬ 
formance of <p\ over 0Trem was not limited to this problem, 
but was present in all Plummer sphere tests with primordial 
binaries we tested. For example, if we take a 5 body Plum¬ 
mer model as described previously in this section, run for 
t/fcross = 2.9, and calculate 5 = \AE/E\, <p\ yields f cpu = 64 
with 8 = 9.6 x 10~ 6 while (pTiem yields t cpu = 100 and 
8 = 1.3 x 10 -5 : better conservation for smaller compute 
times. 


4 SYMPLECTIC AND REVERSIBLE 

INTEGRATOR COMBINING KICKS AND 
KEPLER SOLVERS 

(ph uses Kepler solutions in place of the faster kick steps of 
DKD leapfrog. We now describe a symplectic and reversible 
integrator combining the speed of leapfrog and the collisional 
accuracy of Kepler solvers. To do this, we write a first order 
non-reversible symplectic integrator as we did in eq. 01. Let 
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Figure 4. Accuracy achieved for a given computational effort. 5 particles are sampled from a Plummer sphere and a binary is included 
as described in the text. The binary is twice as bound as the Plummer sphere in standard units. Numerical results are shown for the 10 
integrals of motion including 6 center of mass phase space variables, outperforms the other methods. 


A be a set of pairs of particles and A c is its complement 
(all the other pairs). Then the pseudocode for <f>' h is shown 
in Algorithm [3j With operators, we have eq. (|15[) : 

<j>h = fl ex P (hDK i: ) exp (-hDr tj ) I exp (hD w ) exp (/ 
\(,id)eA ) 

Yli’h’i’-lj i’U’l, 

(15) 

with 

(i,j)eA c 


Drift all particles for time h ; 

for pairs of particles (i, j) in A c do 

Kick particles i and j using their mutual force 
only, 
end 

•t) • for pairs of particles (i, j) in A do 

Drift particles i and j for time —h ; 

Apply a Kepler solver to advance the relative 
coordinates of i and j by h ; 

Advance center of mass coordinates of i and j by 
h ; 
end 

Algorithm 3: Pseudocode for </>' h 
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Figure 5. Energy conservation for a given computational effort 
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F igur e 6. applied to the f igure-eight three body problem 
(jChenciner & Montgomervll2000h . A time step h = 0.001 is used 
and the problem is run for 100 periods. The energy error is aver¬ 
aged over the period as in Figure [3] Consistent with existence of 
a nearby H , and periodic behavior characteristic of a symplectic 
integrator, there is no secular drift in the energy error. 


With method cf>' h we could choose, for example, to kick pairs 
that are far from each other and apply the Kepler solver to 
tight binaries. In practice tight binaries will interact with 
other particles and perhaps break up, and this integrator 
may not be practical for such situations. A second order 
method is developed following the step of eq. ©: 

<t>k = 0h/20h/2- (16) 

As we showed in Section a kick step is symplectic and 
conserves 9 of 10 integrals of motion to machine precision. 

We first show the long-term conservation properties 
characteristic of symplectic and reversible methods of <f>%. 
We apply the integrator to the figure-eight three body prob¬ 
lem. We assign the three particles to numbers 1-3. We let 
A = {(1, 2), (1, 3)}, where (1,2) and (1,3) are the pairs of 
particles evolved with Kepler solvers, and A c = {(2, 3)}. We 
choose h = 0.001 and run the problem for 100 periods. The 
result is shown in Figure [5] and we indeed see the excellent 
conservation properties we expect from symplectic and re¬ 
versible methods. We have again averaged the energy error 
over each period. 


Now we apply the integrator to a more complicated 
problem. We take the same figure-eight problem, but now 
replace each particle by a binary whose center of mass be¬ 
haves identically as in the original case. The binary pairs 
will have e = 0.9 and a = 0.01 again as in the case of Figure 
[4] The binaries have total energy 29 times greater in magni¬ 
tude than the energy of the original problem. They complete 
more than 1000 orbits each per each period of the original 
problem. In this case it is not clear whether the KAM theo¬ 
rem will apply, but numerical experiments indicate that for 
a modest number of periods, the system remains close to 
periodic. 

We measure the computational effort required to 
achieve a given accuracy in the integrals of motion. The sys¬ 
tem is run for 2 periods. The result is shown in Figure [TJ In 
this case (j>h performs slightly better than in energy con¬ 
servation, indicating it is also a good collisional integration 
tool. Conservation of other integrals is at machine precision, 
as is the case of <)>^; we can show this through similar analysis 
to that of map K in section m 


5 CONCLUSION 

The main goal of this paper is to present a new colli¬ 
sional A-body symplectic method and its various proper¬ 
ties and perform tests on toy problems. In the small prob¬ 
lems we present, it performs equally as well or better than 
the standard Herm i te me thod and the symplectic method o f 
IPreto fe Tremainel d 1999l l and iMikkoIa fc Tanikawal (il999il . 
We discuss, prove, and numerically show many of the inte¬ 
grator’s properties: its symplecticity, its exact conservation 
of 9 out of 10 integrals of motion, its reversibility, its order of 
integration, and its suitability for collisional problems. We 
perform tests on regular three body problems and chaotic A- 
body problems. Although we have only presented a second 
or der integ r ator, the order can be increased by the method 
of lYoshidal d 199(T) . 

A secondary goal is to present a method such that not 
every pair needs to be treated via Kepler solver, and as a 
result we can increase speed for a given accuracy in some 
problems. 

The next step, which we will take in forthcoming work, 
will be to carry out larger A tests on more realistic prob¬ 
lems. We are especially optimistic about the tests because of 
results bv lGoncalves Ferrari et all d2014l l. since all our tests 
show similar or better performance than SAKURA. The in¬ 
tegrator should prove a useful tool in the study of collisional 
A-body dynamics. 
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APPENDIX A: TIME REVERSIBILITY OF SAKURA 

Here we show that, for N > 2, the SAKURA algorithm presented in Listing 1 of lGoncalves Ferrari et aD (12014! ') is not time 
reversible. We then provide an explanation for the irreversibility. 

Initial conditions at time t are y = {ri,Vi}. r\; and v-, are position and velocity vectors, respectively, for particle i. One 
timestep transforms these initial conditions into final conditions at time t + h: 

r'i = n + hn + y— [K(n — r k , in — v k ,mi + m k , h) - (n + hm) + ( r k + hv k )] , 
k^i 

v'i = Vi+’Y' —- \K(n - r k ,Vi - Vk,rrii + m k ,h) - Vi + v k \ . (Al) 

“ m i + mk L J 

Equation m corresponds exactly to Listing 1. Here, K(r,v,m,t) = x(t) is the solution to the Kepler problem 


d 2 x 
dt 2 


Gmx 


subject to initial conditions x = r, dx/dt = v at t = 0. The Kepler solution itself is time reversible; the question is whether 
algorithm ED is. 

To determine time reversibility, form the pair differences: 


K{n — rj,Vi - Vj,rrn + rrij, h ) 

+ T -—- [K(n - r k , Vi - v k , mi + m k , h) - (n + hn) + (r k + hv k )] 

, mi + m k 

_ - ink - [K(rj - r k , Vj _ v k ,mj +m k ,h) - (rj + hvj) + (r k + hv k )] , 

. mj + m k 

K (n — rj,Vi — Vj,mi + mj,h) 

+ T ——— \K(n - r k ,Vi - v k ,mi + m k , h) - Vi + vA 

. mt + m k L J 

- V - — - \ K(rj - r k ,Vj - v k ,mj + m k ,h) - Vj + v k \ ■ 

. m, + m k L J 


(A2) 


These difference vectors are redundant coordinates and do not form a complete coordinate system. We can complete the 
coordinate system using the center of mass position and the total momentum. We can show using Equations eu that 
SAKURA conserves exactly the center of mass integrals. This result is supported numerically by Figures [2] and [4] It follows 
SAKURA’s center of mass motion is reversible. Consider the case N = 2, for which 

r[ - r' 2 = K(n - r 2 ,vi - v 2 ,mi + m 2 , h) , 

v[ — v 2 = K(m — V2, vi — V2, mi + m2, h) . (A3) 

Clearly, this is time reversible, 

ri — r 2 = K(r[ - r 2 ,v[ - v 2 ,mi+m 2 ,—h) , 

Vx-v 2 = k(r' 1 -r' 2 ,v[-v 2 ,m 1 + m 2 ,-h) . (A4) 


Thus the total solution is reversible. It is also symplectic because it solves the two-body problem exactly. 
However, N > 2 is not time reversible, for example, N = 3: 


r i - r 2 


/ / 
r>i - v 2 


K(r i — r2, vi — V2, mi + m 2, h) 
m 3 

mi + m3 
m 3 

m 2 + m3 
K(r 1 - r 2 ,v 1 — v 2 ,mi +m 2 ,h) 
m3 


[K(m - r 3 ,v 1 - v 3 ,mi + m3, h) — {ri + hvi) + ( r 3 + hv 3 )] 
[K(r 2 - r 3 ,v 2 - v 3 ,m 2 + m3, h) — ( r 2 + hv 2 ) + ( r 3 + hv 3 )] 


+ 


mi + m3 
m3 


m 2 + m3 

While this is more complicated than 


K(m — r 3 , vi — V3 , mi + m3, h) — vi + vaj 
\k(r 2 - r 3 , v 2 - v 3 , m 2 + m. 3 , h) - v 2 + v 3 j . 


(A5) 


it is still possible that the algorithm is time reversible; to be sure, we must evaluate 
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the other two pairs and then explicitly test time reversal. 
/ / 

r l ~ r 3 = 

+ 


K(r i - r 3 ,v i - i) 3 ,mi + m 3 ,h) 
m 2 


mi + m2 
m 2 


[K{ri — r 2 ,vi — v 2 , mi + m2, h) — (n + hv 1) + (r 2 + /1W2)] 
[JC(r 2 - r 3 , v 2 - v 3 ,m 2 + m3, h) — (r 2 + hv 2 ) + ( r 3 + hv 3)] 


Vi-V 3 = 


+ m 2 + m 3 

iT(ri — r 3 ,wi - « 3 ,mi +m 3 ,h) 
m 2 


mi + m 2 
m 2 


[k{r 1 — r 2 , vi — v 2 , mi + m 2 , ft) — ui + V 2 J 
K(r 2 — r 3 , D 2 — i> 3 , m 2 + m 3 , h) — v 2 + v 3 \ 


r 2 -r 3 = 


v 2 - v 3 


m 2 + m3 
i>f(r 2 - r 3 ,v 2 - v 3 ,m 2 + m 3 ,h) 
mi 


+ - 


m 1 + m 2 
mi 


[K(n - r 2 ,v 1 - V 2 ,mi + m 2 , h) — (n + hv i) + (r 2 + hv 2 )\ 
[K(r ! - r 3 ,v 1 - v 3 ,mi + m3, h) — (ri + hv 3 ) + ( r 3 + hv 3)] 


mi + m3 
K(r 2 - r 3 ,v 2 - v 3 ,m 2 +m 3 ,h) 
mi 


+ 


mi + m 2 
mi 

mi + m 3 


K{n — r 2 , vi — v 2 , mi + m 2 , h) — vi + « 2 j 
[K(ri — t- 3 , vi — v 3 ,mi + m 3 , h) - vi + r 3 ] ■ 


(A 6 ) 


The presence of the additional pair terms on the right-hand sides of equations (lA5l) and (1A61) breaks time reversibility. 

The reason why SAKURA is not time reversible is that the Kepler steps in ED are not done in serial but instead are 
done in parallel, using the initial phase space coordinates for every Kepler evaluation. When N > 2, so that more than one 
pair needs to be advanced, the phase space variables need to be updated after each Kepler evaluation. This is a necessary 
but not sufficient condition for time reversibility. An additional requirement is that the pairs be evaluated in a reversible 
order, which is accomplished in our paper by splitting the timestep in half and using one pair order for the first half, and the 
reversed pair order for the second half (this is what is meant by the adjoint map). 
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